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ABSTRACT 


This  thesis  analyses  data  of  and  builds  a  simulation 
model  for  the  track  of  an  underwater  vehicle  as  perceived 
by  a  test  range  of  three  dimensional  short  baseline  sonar 
arrays.  In  this  way  many  random  replications  of  track 
become  available  quickly  and  inexpensively.  These 
simulations  support  a  larger  project  whose  object  is  to 
monitor  the  performance  of  the  test  range  and  provide  clues 
for  troubleshooting  problems.  In  particular,  joint  values 
of  sensor  array  estimated  displacement  and  reorientation 
corrections  are  generated  and  their  error  distribution  is 
quantified.   ^ 
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I.  INTRODUCTION 

The  Naval  Undersea  Meapons  Engineering  Station  (NUWES) 
operates  nine  test  ranges  that  serve  a  variety  of  purposes 
•for  testing  and  validating  the  Navy's  current  and  future 
weapons  systems.  This  thesis  is  in  support  of  a  larger 
project  whose  goal  is  to  monitor  the  performance  of  the 
short  baseline  ranges  (ie.  those  ranges  in  which  each  array 
produces  a  three  dimensional  track  of  moving  vehicles). 
More  specifically,  it  deals  with  the  development  of  a 
computer  simulation  which  generates  realistic  random 
replications  of  an  underwater  vehicle's  track  (ie.  vehicle 
position  perceived  by  a  sensor  array  on  the  range  at  a 
number  of  equally  spaced  time  points).  Such  simulations 
serve  to  provide  information  about  the  inherent  variability 
in  the  output  of  the  tracking  range,  especially  for 
assessing  the  calibration  error.  Its  use  can  provide  an 
improved  understanding  of  the  processes  involved  and  a 
reduction  in  the  frequency  of  range  shutdown  for  the  purpose 
of  resurveying  the  remote  sensor  locations. 

When  the  simulated  replicate  tracks  a^re  passed  to  the 
program  KEYMAIN  (a  project  FORTRAN  program  that  estimates 
displacement  and  orientation  corrections  to  the  range  remote 
sensor  arrays),  the  output  is  used  to  generate  bivariate 
scatter  plots  of  these  corrections.   Although  extensive  work 
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of  this  type  was  not  possible  in  the  current  thesis,  the 
limited  results  indicate  a  surprising  amount  o-f  variability 
and  suggest  that  the  nature  and  extent  o-f  variability  can 
change  noticeably  with  relatively  minor  changes  in  the 
localized  conditions.  Considerable  testing  using  this 
simulation  tool  is  clearly  indicated. 

The  research  reported  here  involves  three  general 
activities: 

(a)  Analysis  o-f  real  underwater  track  data  to  learn  their 
important  behavioral  characteristics. 

(b)  Simulation  model  -formulation  and  construction. 

<c)  Development  and  programming  of  algorithms  to  merge 
with  the  existing  project  programs  to  produce  any 
number  of  random  replications  with  correction 
estimates  for  each  simulated  track. 

The  organization  of  the  thesis  is  as  follows: 
Section  II  contains  explicit  background  material  and 
provides  a  framework  for  the  research  and  shows  how  it 
relates  to  the  larger  project.  Section  III  explains  the 
data  analysis  procedures  and  results.  The  simulation  model 
is  developed  in  Section  IV  and  the  interfacing  with  the 
overall  project  programs  is  described  in  Section  V.  Results 
and  conclusions  are  detailed  in  Section  VI,  with  areas  for 
future  work  listed  in  Section  VII. 

A  number  of  appendices  are  included  to  provide  the 
detailed  support  too  extensive  to  be  included  in  the  main 
body  of  the  paper.  They  are  referenced  in  the  appropriate 
sections.     Additionally,    to   speed   the   completion   and 


availability  of  the  KEYMAIN  program,  the  author  wrote  two 
subroutines  (CONECT  and  REDUCE)  to  be  included  in  that 
package.  Appendices  E  through  G  contain  the  development  o-f 
these  subroutines  and  the  FORTRAN  code  listings. 
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II.  BACKGROUND 

Consider  a  three  dimensional  underwater  tracking  array 
composed  of  four  hydrophones  arranged  to  define  a  local 
Cartesian  coordinate  system  (see  Figure  1).  Such  an  array 
is  capable  of  tracking  a  target  downrange,  crossrange  and  in 
depth  in  its  local  coordinate  system;  it  is  termed  a  short 
baseline  array  because  of  the  short  distance  of  the  X,  Y, 
and  Z  hydrophones  from  the  "corner"  hydrophone  (30  feet) 
that  is  used  to  gather  the  three  dimensional  tracking 
information. 

The  arrays  take  "fixes"  of  target  position.  A  fix  is 
defined  by  a  bearing  (azimuth  and  elevation)  and  distance 
from  the  array  to  the  target.  Tracked  targets  on  the  range 
a^re  fitted  with  an  acoustical  device,  called  a  "pinger", 
that  emits  pulses  of  sound,  or  "pings",  at  a  specific 
frequency  at  precisely  timed,  regular  intervals.  The 
elapsed  time  from  the  generation  of  the  ping  at  the  target 
to  the  reception  of  the  ping  at  the  array  establishes  a 
distance   from   the   target  to  the  Array.  Because   of   the 

distance  that  separates  the  individual  hydrophones  of  the 
array,  the  ping  arrives  at  each  hydrophone  at  a  slightly 
different  time.  This  time  difference  can  be  resolved  into  a 
bearing  from  the  array  to  target.  (Actually,  these  fi>;e5 
are  "apparent"  fixes.    They  are  used  to  initialize  a   sound 
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Figure  1  :  3-D  Sensor  Array 
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ray  tracing  algorithm  that  leads  to  the  "real"  fix.)  A 
series  o-f  such  fixes  establishes  a  target  track  in  that 
particular  array's  own  local  Cartesian  coordinate  system. 
If  the  precise  position  of  the  array  is  known,  its  local 
track  information  can  be  translated  to  one  common  range 
coordinate  system. 

Position  of  an  individual  array  in  the  range  Cartesian 
coordinate  system  is  described  not  only  in  terms  of  its 
downrange,  crossrange  and  depth  (X,Y  and  Z)  coordinates 
(termed  location),  but  also  with  respect  to  X-tilt,  Y-tilt 
and  Z-rotation  (commonly  called  roll,  pitch  and  yaw)  angles 
from  the  range  coordinate  system  axes  (termed  orientation). 
Both  sets  of  measures  are  needed  to  translate  accurately 
from  local  to  range  coordinates. 

Typically,  a  range  is  composed  of  many  individual 
arrays.  Nanoose  Range,  for  example,  has  24  arrays,  while 
Dabob  Bay  has  7.  The  arrays  are  arranged  in  such  a  way  that 
the  array  coverages  overlap  one  another  to  provide 
continuous  tracking  on  a  target  vehicle.  Such  overlapping 
areas  a^re  called  crossover  regions  (see  Figure  2)  ,  and 
produce  two  sets  of  track  on  the  same  target  for  the  same 
time  period. 

Ideally,  the  corresponding  points,  or  fixes  described 
earlier,  from  each  array  in  a  crossover  region  should 
translate  to  identical  points  on  the  range  coordinate 
system.   In  practice,  this  seldom  happens. 
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Figure  2  :  Crossover  Regions 
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Three   major   sources   of  this  variability  in   crossover 

track  data  are: 

(a)  Slippage   o-f   the   sensor  arrays  -from    their   assumed 
positions  in  the  range  coordinate  system. 

(b>  Time  synchronization  and  instrumentation  problems. 

(c)  Inhomogenei  ties   and    temporal   variability   in    thie 
water  column. 

Although  the  reasons  -for  slippage  of  the  arrays  are 
speculative,  the  fact  that  slippage  occurs  is  evidenced  by 
the  change  in  sensor  locations  after  a  range  resurvey  is 
performed.  The  question  of  discriminating  timing  errors, 
item  (b) ,  from  slippage  errors  is  on  a  future  agenda. 
Investigation  done  by  Main  CRef.  13  provides  us  with 
methodology  for  treating  the  random  components  of  these 
errors.  Item  <c)  falls  into  a  broader  category  which  may 
require  extensive  investigation.  It  may  also  provoke  a 
reassessment  of  the  assumed  uniform  horizontal  quality  of 
the  range  water  column.  It  seems  wise  to  account  for 
slippage  and  instrumentation  problems  first.  Our  work  is 
confined  to  (a). 

The  present  research  is  in  direct  support  of  Professor 
Robert  R.  Read  of  the  Naval  Postgraduate  School  who  has  been 
working  on  the  slippage  question.  He  has  produced  a  FORTRAN 
program  (KEYMAIN)  which  takes  the  crossover  data  and 
estimates  array  position  corrections  using  a  non  linear 
least  squares  algorithm.  The  surfaces  that  enter  into  this 
least  squares  optimization  function  are  rather  flat  and   the 
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values  do  not  change  much  as  the  location  corrections  are 
varied.  Replicated  data  is  needed  to  quanti-Fy  the  extent  to 
Mhich  the  estimated  corrections  are  Mi  thin  the  range  of 
natural  variability.  The  simulation  model  in  this  thesis 
provides  the  tool  by  Mhich  this  variability  can  be 
quantified.  Mith  the  simulation,  the  needed  replicated  data 
can  be  quickly  and  easily  generated,  and  scatterplots  of 
displacement  (change  in  location)  versus  rotation  (change  in 
orientation)  can  be  made  to  investigate  the  inherent 
variability  in  the  correction  estimates. 
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III.  DATA  ANALYSIS 

A.   GENERATING  RESIDUALS  FROM  REAL  TRACK  DATA 

The   real  track  data  used   for  this  study  was  taken  -from 

Nanoose   range   in   September,   1982   (Figure   3).    It   was 

supplied   as  an  N  X  7  matrix,   N  representing  the  number   of 

data   points   in   the  crossover  data   set   and   the   columns 

representing  the  -following: 

Column  1  Point  count.  Every  pulse  emitted  by  the 
pinger  is  assigned  a  sequential  number 
beginning  at  the  start  o-f  the  tracking  run. 
Therefore,  this  column's  value  increases  by 
one  each  row  unless  a  data  point  is  missing, 
which  would  be  indicated  by  a  sequential 
omission  in  this  -first  column. 

Columns  2-4  The  X,Y,  and  Z  coordinates  o-f  the  target 
■for  the  column  one  point  count  as  determined 
by  the  first  sensor  in  the  crossover  pair. 

Columns  5-7  The  X,Y,  and  Z  coordinates  of  the  target  for 
the  column  one  point  count  as  determined  by 
the  second  sensor  in  the  crossover  pair. 

There  were  six  such  data  sets  used.  These  particular 
sets  were  chosen  because  each  formed  a  straight  line  in  the 
range  space  and  so  the  best  straight  line  path  of  the  target 
could  be  estimated  from  the  data.  (In  contrast,  estimating 
the  best  curved  path  from  curved  track  data  would  have  been 
far  more  challenging,  but  with  little  or  no  gain  in  the 
examination  of  residuals.) 

The  idea  was  to  take  the  track  data,  fit  the  best 
straight   line  possible  through  the  data,   then  examine   the 


17 


CM 
00 

CL 
LlJ 
CO 

CD 
I 


LU 
O 

< 

UJ 

cn 
O 
O 


o 
< 


O  < 

^  q::  CD  + 


(/) 


< 


0002 


o 

Ld 
if) 


g 


< 

Cr   IT)    + 

< 


^  s  ^ 

LU    << 

ic  n    ^ 


5 
o 

to 


+ 


o 

o 

o 

<r) 

lO 

cr 

•  H 

o 

cr 

o 

-P 

o 

(— 

cr. 

lO 

UJ 

LJ 

0) 

Li_ 

tL 

^"^ 

C 

LU 

c- 

O 

pr. 

z 

< 

G.' 

Q^ 

c; 

-z. 

C 

5 

c 

o 

o 

C 

o 
o 

m 

Q 

o 
o 
o 


0003- 


(133J)  30NVdSS0dD 


18 


residuals  formed  by  subtracting  the  fitted  straight  line 
point  sequences  from  the  original  data. 

The  commented  FORTRAN  program  written  to  generate  the 
residuals  is  included  as  Appendix  A,  and  its  mechanics  are 
discussed  below. 

Since   the   data   was    given  in   a   single   matrix   and 

residuals   were  desired  for  each  sensor,   the  data  was  first 

separated  into  two  tracks,   each  an  N  x  3  matrix;   call  them 

TRl  and  TR2.   Next,  a  3  x  3  covariance  matrix  for  each  tracit: 

was  computed  by  the  formula 

N  

GOV  =  C  ZH  (TR(i)-TR)'  (TR<i)-TR)3  /  N-1 
i  =  l 

where 

TR(i)   =3  component  row  vector  of  the  X,Y,Z  coordinates  of 

the  i    data  point 

TR  =  3   component   row  vector  of  the  average  of   the   N 

<X,Y,Z)  coordinates  for  that  track  data 

N   =   number  of  data  points 

The   best  straight  line  estimate  of  the   target  path   is 

the  straight  line  such  that  the  sum  of  the  distances  of  each 

data  point  to  the  line  is  a  minimum.    This  implies  that  the 

distance  from  each  point  to  the  line  is  a  minimum  -  that  is, 

each   point  should  be  projected  onto  the  line   orthogonally. 

The  method  employed  to  accomplish  this  orthogonal  regression 

was   that   of   principle   components.   Principle   components 

requires   that   the  eigenvector  associated  with  the   largest 
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eigenvalue  o-F  the  covariance  matrix  be  identi-fied.  It  Mas 
computationally  convenient  to  make  a  3  x  3  matrix  o-f  the 
eigenvectors,  Mith  the  -First  column  being  the  eigenvector 
associated  Mith  the  largest  eigenvalue,  the  second  column, 
the  eigenvector  associated  Mith  the  second  largest 
eigenvalue  and  third  column,  the  eigenvector  associated  with 
the  smallest  eigenvalue.  This  done,  the  -folloMing  matrix 
mul t  i  p 1 i  cat  ion 


PROJECTION  =  (TR  -  AVE)  x  EIGENVECTORS 
where 

TR  =  N  X  3  track  data  matrix 


AVE  =  N   X   3  matrix  made  o-f  N  identical  rows   o-f 
the   X,Y,  and  Z  averages  o-f  the  track  data 
EIGENVECTORS   =3X3  eigenvector  matrix  described  above 
yields  the  N  x  3  matrix  PROJECTION,   whose  entries   are   the 
orthogonal   projections   o-f  the  track  data  points   onto   the 
axes   o-f   a  new  3  dimensional   Cartesian   coordinate   system 
whose   origin  is  located  at  TR  (the  (X,Y,Z)  averages  -for  the 
data   set)   and  whose  axes  are  rotated  such  that  the   new   X 
axis   is  the  best  straight  line  that  describes  target   path. 
For   example,   the   -first   row   o-f   PROJECTION   is   a   three 
component   row  vector;   call  the  components  (pl,p2,p3).  Then 
the  point  (pi, 0,0)  is  the  projection  o-f  the  -first  track  data 
point   onto   the  new  X  axis  (known  to  be  the   best   straight 
line  that  describes  the  target  path),   the  point  (0,p2,0)  is 
the  projection  o-f  the  -first  track  data  point  onto  the  new   Y 
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axis,  and  the  point  (0,0, p3)  is  the  projection  of  the  first 
track  data  point  onto  the  new  Z  axis.  Since  we  Are  only 
interested  in  the  projections  onto  the  new  X  axis,  we  can 
replace  the  second  and  third  columns  o-f  PROJECTION  with 
zeros  and  have  the  coordinates  of  the  best  straight  line 
path  of  the  target.  These  coordinates  are  in  the  PROJECTION 
coordinate  system,  and  must  be  translated  back  into  the 
range  coordinate  system.  Using  the  matrix  PROJECTION  (with 
last  two  columns  =  0)  and  performing  the  fallowing 
mul tipl i cat ion 

OLD  =  EIGENVECTORS  x  PROJECTION' 
yields  the  3  x  N  matrix  of  the  track  data  points  in  the 
range  coordinate  system.  Note  that  OLD  must  be  transposed 
to  get  the  N  x  3  format  desired.  Since  the  mean  had  been 
subtracted  off  before  computing  the  PROJECTION  matrix,  to 
get  the  data  points  back  to  the  proper  position  requires 
that  the  means  be  added  back  in.   The  matrix  addition 


OLD'  +  AVE 
yields  the  N  x  3  matrix  of  the  orthogonal  projections  of  the 
data   onto   the   straight   line  path  of  the   target   in   the 
original  range  coordinate  system. 

Actually,  after  having  computed  the  PROJECTION  matrix, 
there  is  yet  another  step  to  take  before  translating  back 
into  the  range  coordinate  system.  Based  on  an  assumption  of 
constant  speed  for  the  target,  the  track  data  points  should 
be   equally   spaced  since  the  pulses  are    emitted  at  regular, 
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precisely  timed  intervals.  This  could  be  accomplished  by 
using  an  interval  equal  to  the  total  distance  divided  by  the 
number  o-f  data  paints,  but  i-F  data  points  are  missing  in  the 
track  data  set  (and  each  set  o-f  track  data  Mas  missing 
several  points) ,  then  this  method  introduces  an  error.  The 
method  chosen  should  shoM  that,  i-f  a  single  data  point  is 
missing,  the  gap  between  the  consecutive  points  on  the  data 
track  should  be  twice  as  long  as  i-F  none  were  missing.  The 
■first  column  (  the  point  counts)  o-f  the  N  x  7  track  data 
matrix  contains  the  in-formation  on  missing  points. 
Performing  a  simple  least  squares  regression  of  the  point 
counts  onto  the  first  column  of  the  PROJECTION 
matrixtranslates  this  information  into  the  best  straight 
line  point  sequence  for  target  motion.  Now  when  the 
projections  onto  the  first  principle  component  are 
translated  back  onto  the  range  coordinate  system,  the 
result  is  truly  the  best  set  of  track  points  attainable 
using  all  the  information  contained  in  the  track  data  set. 

To  obtain  the  residuals,  the  straight  line  estimated 
track  points  are  subtracted,  point  for  point,  from  the 
original  track  data.  It  was  important  to  discover  the 
distribution  of  the  residuals  for  each  track  segment  because 
to  simulate  the  track  segments  later,  residuals  were 
simulated  and  then  added  to  the  straight  line.  Knowing  the 
character  of  the  residuals  allowed  the  generation  of  very 
realistic  track  data  for  the  simulation  model. 
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B.  ANALYZING  TRACK  DATA  RESIDUALS 

Each  o-f  the  six  track  segments  produced  two  sets  of 
residuals  —  one  set  -for  each  o-f  the  two  sensors  of  the 
crossover  pair.  Each  set  was  further  divided  into  X,  Y,  and 
Z  components.  There  were,  therefore,  36  sets  of  residuals  to 
be  examined  for  distributional  characteristics.  Graphical 
analysis  was  performed  on  each  set  of  residuals  (histograms, 
cumulative  density  plots  and  QQ  plots)  and,  when,  from  that 
analysis,  a  distribution  for  the  residuals  could  be 
determined,  formal  statistical  tests  were  performed  to 
verify  that  the  best  distribution  was  chosen.  The  analysis 
showed  the  residuals  to  be  normally  distributed.  The  best 
and  worst  case  graphical  and  analytical  results  are  included 
in  Appendix  B.  Note  that  there  were  some  very  good  fits  to 
a  normal  density  (pp.  66-71)  with  statistical  significances 
for  the  Chi  Squared  and  Kolmogorov-  Smirnov  goodness  of  fit 
tests  well  above  a  very  conservative  .35  in  all  but  one 
case.  Even  in  the  data  that  most  poorly  resembled  a  normal 
density  (pp.  72-77) ,  the  Kolmogorov  -  Smirnov  goodness  of 
fit  significance  level  never  fell  below  .25.  From  tins 
analysis  it  was  concluded  that  a  normal  density  for  the 
residuals  was  accurate  and  appropriate  for  the  simulation 
model . 

C.  TIME  SERIES  ANALYSIS 

From  a  simulation  point  of  view,  it  is  very  convenient 
to   assume  independence  of  the  residuals  between   successive 
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points   in   the  original  track.    I-f  the  assumption   is   not 

made,   then   one  must  resort  either  to  using  a   point   count 

interval   o-f  some  integer  value  greater  than  one  for   making 

the    simulated   track   -from   the   original,    or   build   an 

autoregressi ve  model. 

Concern  -for  this  correlation  o-f  the  data  over   time   led 

to  a  time  series  analysis  o-f  the  residuals.     This  was  done 

in   two   ways  -for  each  track  -from  a  single   sensor.    Recall 

that  for  each  data  point,  residuals  were  generated  in  the  X, 

Y,  and  Z  directions  o-f  the  range  coordinate  system.   Each  o-f 

these  components  was  subjected  to  a  time  series  analysis   to 

determine   whether   there  was  a  dependence  in  the  errors   in 

any  single  direction  over  time.    Then,   the  distance  o-f  the 

data  point  -from  the  straight  line  target  path,  given  by 

SQRTC  (RES  ^)  +  (RES  ^)  +  (RES  ^)  : 
X  y  z 

was  examined  to  determine  i-f  the  magnitude  o-f  the  error  of 
one  point  was  correlated  with  the  magnitude  of  the  error  of 
the  next  point.  The  time  series  analysis  examines  the 
dependence  from  point  to  point,  on  every  other  point,  on 
every  third  point,  and  so  on,  so  that  the  interval  at  which 
one  can  assume  independence  (indicated  by  an  autocovariance 
value  of  zero)  can  be  determined.  Appendix  C  contains  the 
autocorrelation  graphs  for  the  6  data  sets;  first  the  error 
magnitudes  are  examined  for  both  sensors  of  a  data  set, 
followed  by  the  individual  error  analysis  in  the  X,  Y,  and  Z 
directions. 


24 


There  Mas  no  clear  cut  interval  -for  which  independence 
seemed  to  emerge  in  the  data  sets.  There  is  instead  a 
random  pattern  o-f  insignificant  dependence  irom  the 
beginning  o-f  the  autocorrelation  analysis,  leading  to  the 
conclusion  that  point  to  point  independence  was  appropriate 
•for  the  simulation  model.  Attention  is  directed  to  the 
noise  in  the  autocorrelation  graphs  with  no  distinct 
patterns  emerging,  and  correlation  magnitudes  smaller  than 
0.2  in  most  cases.  It  was  decided,  there-fore,  that  any 
important  time  correlation  was  not  so  large  as  to  cause 
excessive  or  even  noticeable  error  in  the  simulated  tracks. 

D.   STUDY  OF  RESIDUALS 

The  study  o-f  the  residuals  yielded  important  interesting 
in-f  ormation ,  summarized  in  Table  1.  The  -first  column  oi 
Table  1  is  the  standard  deviation  o-f  the  residuals  from  the 
le-ft  sensor  in  the  downrange  (X),  crossrange  (Y)  and  depth 
(Z)  directions,  three  values  -for  each  data  set.  Column  8 
gives  the  correspondjing  in-formation  -for  the  right  sensor  o-f 
the  crossover  data  pair.  Columns  2,  3  and  4  -for  each  data 
set  are  the  columns  o-f  a  3  x  3  correlation  matrix  for  the 
residuals  o-f  the  le-ft  sensor,  columns  5,  6  and  7, 
corresponding  in-formation  for  the  right  sensor. 

First,  note  the  difference  between  the  left  and  right 
sensor  standard  deviations  of  residuals  in  all  three 
directions.  Although  there  are  some  very  good  comparisons 
(data   set  D2A1,   crossrange,   and  data  set  D5A ,   downrange) 
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TABLE  1 
STANDARD  DEVIATIONS  OF  RESIDUALS 
AND 
CORRELATION  COEFFICIENTS  OF  RESIDUALS 


SD 
LEFT 


CORRELATION 
LEFT 


CORRELATION 
RIGHT 


SD 
RIGHT 


DATA  SET  D2A1 


N  =  67 


0.842 

1.000 

0.916   -0.899 

1.000 

-0.863 

0.737 

0 ,  72 1 

2.  141 

0.916 

1.000   -0.728 

-0.863 

1 .  000 

-0.455 

2.  154 

2.394 

-0.899 

-0.728    1.000 
DATA  SET 

0.787 
D2A2  ,  N 

-0.455 
=  70 

1 .  000 

2.063 

1.484 

1 .  000 

0.859   -0.787 

1 .  000 

-0.360 

0.040 

1.  163 

2.396 

0.859 

1.000   -0.759 

-0.360 

1 .  000 

-0.547 

2.  147 

2.694 

-0.787 

-0.759    1.000 

0.040 

-0.547 

1 .  000 

1.976 

DATA  SET  D2B 


N  = 


0 .  693 
2.000 
2.353 


1.000    0.738   -0.888    1.000   -0.949    0.902    0.937 

0.783    1.000   -0.434   -0.949    1.000   -0.728    3.350 

-0.388   -0.434    1.000    0.902   -0.728    1.000    3.042 


DATA  SET  D4 


N  =  9Z 


0.497 

1 .  00.0 

0.741 

-0.904 

1 .  000 

-0.681 

0.380 

0.  502 

1.955 

0.741 

1 .  000 

-0.687 

-0.631 

1 .  000 

-0. 146 

2.467 

1.537 

-0 . 904 

-0.687 

1 .  000 

0 .  380 

-0. 146 

1 .  000 

1.843 

DATA  SET  D5A 


N  =  3: 


1.731 

1 .  000 

-0.930 

-0.429 

1 .  000 

-0. 953 

0.163 

5.751 

-0.930 

1  .  000 

0.  152 

-0 . 953 

1 .  000 

-0. 148 

1.338 

-0.429 

0.  152 

1 .  000 

0.  163 

-0. 148 

1 .  000 

1.736 
6.  645 
1.378 


DATA  SET  D5B 


N  =  37 


0.386 

1 .  000 

-0.408 

-0.450 

1 .  000 

0.445 

0 .  395 

0.717 

2.931 

-0.403 

1 .  000 

-0.454 

0.445 

1 .  000 

-0.279 

1.613 

3.397 

-0 . 450 

-0.454 

1 .  000 

0.395 

-0.279 

1  .  000 

4.254 
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there  are  also  some  wide  disparities  (data  set  D2B, 
crossrange  and  data  set  D5B,  crossrange) .  Using  the 
variance  ratio  test  described  in  Larson  CRef.  2:  pp.  449],  a 
statistical  test  Mas  performed  to  determine  whether  the 
variances  (the  recorded  standard  deviations,  squared)  of  the 
left  and  right  sensor  arrays  for  downrange,  crossrange  and 
depth  were  the  same  for  a  given  data  set.  The  results  are 
given  in  Table  2.  Using  0.05  as  the  basis  for  rejection  of 
the  null  hypothesis  (H  :  the  two  variances  are  the  same) ,  8 
of  the  18,  or  44%,  of  the  individual  tests  fail.  With  a 
confidence  level  of  .95,  one  would  expect  roughly  1  failure 
in  the  18  tests.  Eight  failures  is  strong  evidence  that  the 
standard  deviation  figures  do  not  match  very  well. 

Recalling  the  test  data  from  Figure  3  (p.  18),  it  is 
seen  that  there  are  three  sensor  arrays  that  contributed  the 
data:  arrays  4,  5  and  6.  In  data  sets  D2A1 ,  D2A2  and  D2B, 
array  4  is  the  left  array  and  array  5  is  the  right  sensor 
array.  For  data  sets  D4,  D5A  and  D5B,  sensor  array  Sis  the 
left  array  and  array  6  is  the  right  array.  One  might  expect 
that  the  standard  deviations  in  any  single  direction  would 
be  the  same  for  a  given  sensor.  This  turned  out  not  to  be 
the  case.  Using  Bartlett's  test  for  the  equality  of  several 
variances  CRef.  3:  pp.  225-227D,  the  hypothesis  that  the 
variances  for  a  given  sensor  in  a  given  direction  are  equal 
was  tested.  The  results  are  given  in  Table  3.  The  first 
two   sections  of  the  table  compare  the  3  standard   deviation 
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TABLE  2    TEST  OF  THE  HYPOTHESIS  THAT  THE 

STANDARD  DEVIATIONS  WITHIN  A  DATA  SET 

ARE  EQUAL 


VARIANCE 

DEGREES 

LEVEL  OF 

• 

RATIO 

OF  FREEDOM 

SIGNIFICANCE 

DATA  SET  D2A1 

DOMNRANGE 

1.363 

66 

.211 

CROSSRANGE 

.988 

66 

.935 

DEPTH 

1.363 

66 

.229 

DATA  SET  D2A2 

DOWNRANGE 

1.627 

69 

.045 

CROSSRANGE 

1.244 

69 

.364 

DEPTH 

1.858 

69 

.011 

DATA  SET  D2B 

DOWNRANGE 

.547 

52 

.031 

CROSSRANGE 

.357 

52 

.0003 

DEPTH 

.598 

52 

.067 

DATA  SET  D4 

DOMNRANGE 

.982 

92 

.884 

CROSSRANGE 

.628 

92   . 

.027 

DEPTH 

.692 

92 

.079 

DATA  SET  D5A 

• 

DOWNRANGE 

.995 

81 

.939 

CROSSRANGE 

.749 

81 

.  195 

DEPTH 

.508 

81 

.003 
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TABLE    2       (continued) 


DATA  SET  D5B 

DOMNRANGE 

1.527 

86 

CROSSRANGE 

3.300 

86 

DEPTH 

.638 

86 

.051 
0.0 
.038 
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Table  3   TEST  OF  THE  HYPOTHESIS  THAT  THE 
STANDARD  DEVIATIONS  OF  RESIDUALS  OF 
A  PARTICULAR  SENSOR  ARE  EQUAL 

Chi  Squared  Random  Variables 

DEGREES  OF  FREEDOM    .99  QUANT I LE    .999  QUANT I LE 

2  9.21  13.82 

5  15.09  20.52 


DOWN 

CROSS 

DEGREES 

RANGE 

RANGE 

DEPTH 

OF 

FREEDOM 

SENSOR 

4 

(LEFT) 

39.57 

2.02 

1.41 

2 

SENSOR 

5 

(RIGHT) 

14.75 

16.25 

13.96 

2 

SENSOR 

5 

(LEFT) 

131. 14 

101.33 

91.40 

2 

SENSOR 

6 

(RIGHT) 

148-72 

175. 18 

84.31 

2 

SENSOR  5  (LEFT  AND 
RIGHT) 


153.04 


150.87    105.51 


LEFT  SENSORS 
RIGHT  SENSORS 


171.45     162.29     94.48 
168.37     236.97    107.15 


5 
5 
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values  -for  each  sensor  in  a  single  direction.  The  test 
statistic  is  distributed  as  a  Chi  Squared  random  variable 
with  2  degrees  o-f  -freedom.  The  .99  and  .999  quant iles  ior 
such  a  random  variable  are  9.21  and  13. Q2.  In  only  two 
cases  o-f  the  12  trials  would  the  3  standard  deviations  be 
considered  statistically  the  same. 

The  third  section  o-f  Table  3  is  included  because  sensor 
array  5  was  used  as  both  a  right  and  left  array.  Therefore, 
there  were  actually  six  values  -for  downrange,  crossrange  and 
depth  for  that  particular  array.  Using  the  same  test  as 
before,  the  question  of  whether  the  six  values  were 
statistically  the  same  was  investigated.  The  tabulated 
statistic  for  the  3  cases  is  distributed  as  a  Chi  Squared 
random  variable  with  5  degrees  of  freedom.  The  .99  and  .999 
quantiles  for  such  a  ramdom  variable  are  15.09  and  20.52. 
In  every  case,  the  p-value  of  the  test  is  essentially  equal 
to  0.0.   The  results  are  highly  significant. 

Finally,  the  bottom  of  the  table  investigates  whether 
downrange,  crossrange  and  depth  residual  standard  deviations 
are  the  same  for  left  and  right  sensors  in  general. 
Employing  the  same  test  for  the  6  values  produced  the 
tabulated  statistics.  Since  the  test  statistic  is  again  Chi 
Squared  with  5  degrees  of  freedom,  the  level  of  significance 
in  every  case  is  essentially  0.0. 

A  third  question  of  interest  from  Table  1  is  whether 
the   correlation  coefficients  are  the  same  far  the  left   and 
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right  sensors  in  a  data  set.    This  Mas  investigated  in   the 
■folloMing  manner. 

First,  the  normalizing,  inverse  hyperbolic  tangent 
transformation  CRef.  3:  p.  365D  Mas  applied  to  each  o-f  the 
o-f-f  diagonal  correlation  coe-f f iceients  in  the  correlation 
matrices.  (Note  that  the  matrices  are  symmetric,  so  there 
are  only  three  values  o-f  concern  -for  each  matrix.)  This 
trans-formation  makes  each  o-f  the  values  normal  Mith  mean 

1  1  +  rho 
-  In  

2  1  -  rho 

and  variance 


N  -  3 
Mhere  N  is  the  number  o-f  data  points  contributing  to  the 
correlation.  Under  the  assumption  that  the  tMO  independent 
sample  correlation  coe-f -ficients  come  -from  the  same 
population,  their  di-f-ference  is  distributed  normally  Mith 
mean  Q  and  variance 

2 


N  -  3 
We  can  therefore  go  to  the  standard  normal  tables  to   obtain 
the  significance  of  the  statistic 


SORT  C2/(N-3) : 
Mhich    tests   the   hypothesis   that   the   tMo    correlation 
coefficients  are  equal.   These  values  are  tabulated  in  Table 
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TABLE  4   TEST  OF  THE  HYPOTHESIS  THAT  THE 
CORRELATION  COEFFICIENTS  WITHIN  A  DATA  SET 

ARE  EQUAL 


DATA     CORRELATION  STANDARDIZED      LEVEL  OF 

SET     COEFFICIENT     Zl  -  Z2       NORMAL       SIGNIFICANCE 


0.0 

0.0 

0.014 

0.0 

0.0 

0.028 

0.0 

0.0 

0.022 

0.0 

0.0 

0.0 

O.  196 

0.0 

0.057 

0.0 

0.0 

O.  149 


D2A1 

X,Y 

D21A 

X,Z 

D2A1 

v,z 

D2A2 

X,Y 

D2A2 

X,Z 

D2A2 

Y,Z 

D2B 

X,Y 

D2B 

X,Z 

D2B 

Y,Z 

D4 

X,Y 

D4 

X,Z 

D4 

Y,Z 

D5A 

X,Y 

D5A 

X,Z 

D5A 

Y,Z 

D5B 

X,Y 

D5B 

X,Z 

D5B 

Y,Z 

2.867 

16.227 

2.534 

-14.333 

-.434 

-2.453 

1.668 

9.654 

1.  103 

-6.384 

-.381 

-2.203 

2.888 

14.438 

2.896 

-14.481 

.459 

2.293 

1.783 

1 1 . 959 

1.894 

-12.706 

-.695 

-4.659 

.206 

1.293 

-.628 

-3.948 

.302 

1.900 

-.911 

-5.904 

-.903 

-5.850 

-.203 

-1.443 

oo 


4.  To  get  a  signi-ficance  o-f  greater  than  .05  in  this  test 
requires  a  value  between  +/-  1.96,  so  only  3  of  the  IS  pairs 
o-f  correlation  coe-f -f icients  are  statistically  equal,  giving 
very  strong  evidence  that  the  correlation  coefficients  of 
the  residuals  as  recorded  by  the  two  sensors  in  a  crossover 
pair  are,  in  general,  not  equal. 

The  practical  significance  of  the  preceding  tests  is 
that  there  appears  to  be  much  local  variation  in  the  range 
and  that  one  single  model  for  simulating  random  replications 
of  underwater  track  is  not  apparent.  Therefore,  each  case 
must  be  simulated  and  studied  separately,  and  the  study  of 
the  error  estimates  distributions  done  on  a  case  by  case 
basis. 


IV.  SIMULATION  MODEL 

The  simulation  model  itsel-f  is  a  logical  extension  o-f 
the  residual  generation  program.  In  -fact,  the  method  used 
to  simulate  a  data  track  Mas  to  simulate  residuals  and  add 
them  to  the  -fitted  straight  line  obtained  -for  that  track, 
rather  than  to  start  -from  scratch  and  simulate  a  totally  new 
track  at  each  iteration.  This  method  was  very  quick  and 
yielded  very  good  simulated  track  segments.  Figure  4  on  the 
•following  page  is  graphical  comparison  o-f  a  typical  original 
track  segment  overlaid  by  a  track  segment  simulated  by  the 
method  stated  above.  The  comparison  demonstrates  the 
realistic  quality  o-f  the  simulated  track. 

The  regression  method  used  to  compute  residuals  o-f  a 
track  segment  also  produced  the  best  straight  line  in  three 
dimensional  space  to  approximate  the  target  path  through  the 
water.  Given  the  straight  line,  the  residuals  themselves, 
and  the  assumption  of  normality  o-f  the  residuals,  it  is 
possible  to  simulate  a  set  o-f  residuals  -from  normal  (0,1) 
deviates  and  add  them  to  the  straight  line  to  obtain  a 
simulated  track. 

In  simulating  residuals,  the  object  is  to  compute  an 
N  X  3  matrix,  using  normal  (0,1)  deviates,  whose  vectors  o+ 
X,  Y,  and  Z  components  are  normally  distributed  with  mean  O 
and   whose   covariance   matrix  is  equal   to   the   covariance 
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matrix  of  residuals.   That  is 

R  =  T  X  X' 
where 

R  =  3  X  N  matrix  o-f  simulated  residuals 

X  =  N  X  3  matrix  of  N<0,1)  deviates 

T  =  3  X  3  trans-formation  matrix  such  that 

T  X  T'  =  covariance  matrix  o-f  residuals. 
Note   that   R  must  be  transposed  to  get  the  desired   N   x   3 
residual  format. 

It   is   easy  to  show   that  any  3x3  matrix  T   will   not 
alter  the  mean  o-f  0. 

ecr:  =  ECT  X  X': 

=  T  X  ECX  '1 
Since   the   expectation  o-f  X  —  consisting  o-f   normal   (0,1) 
deviates  —  is  identically  0, 

ECR]  =  T  X  (0)  =  0 
giving  the  desired  result. 

The   condition   that  T  x  T'  is   equal  to  the   covariance 
matrix  o-f  the  residuals  is  necessary  because 
COVCR:  =  ECR  X  R'3/N   (the  mean  is  0) 
=  ECT  X  X'  X  X  X  T' ]/N 
=  T  X  ECX'  X   X:/N  X  T' 
(since  T  is  a  linear  operator) 

=  T  X  I  X  T' 
(because  the  covariance  matrix  of  N(0,1)  random  variables  is 
the  identity  matrix) 
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=  T  X  T' 
The   problem  now  becomes  one  o-f  finding  a  3  x  3  matrix   such 
that 

C0VCR3  =  T  X  T' 
We  have  the  covariance  matrix  of  the  residuals.   Let  T  be  an 
upper  triangular  matrix.   Then 

^    "CQV^J   C0Vj2   COVJ 
CQV^j   COV22   COV23 
C0V3J   COV32   COV33 
Noting    that    the   covariance   matrix   is   symmetric    and 
performing  the  matrix  multiplication  yields 

T33  =  SORT  (COV33) 
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^22  =  2°'^^  <^°^22  -  "^23^ 
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^11  =  SQRT<COV^^  -  T^/    -    T^/) 
Since  the  matrix  T  is  easily  computed  and  X  can  be  generated 
from  a  random  number  generating  package,   the  residuals,   R, 
are  quickly  and  economically  computed  for  as  many   simulated 
tracks  as  desired. 
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Adding  the  residuals  thus  -formed  to  the  best  straight 
line  target  path  of  the  original  data  yields  a  simulated 
track. 

The  commented  FORTRAN  program  written  to  per-form  the 
simulation  is  included  as  Appendix  A.  It  is  basically  a 
driver  program  that  generates  a  user  speci-fied  number  o-f 
simulated  tracks  and  interfaces  with  KEYMAIN,  which 
estimates  displacement  and  angular  rotation  values.  Because 
some  of  the  user  friendly  attributes  of  KEYMAIN  interfere 
with  the  speedy  generation  of  the  simulation's  required 
output  parameters,  the  package  has  been  altered  somewhat  to 
increase  speed.  The  output  of  the  driver  program  is  two 
data  files.  One  data  file,  ROTATE. DAT,  i s  an  N  x  4  matrix 
that  gives,  in  the  first  column,  the  maximum  angle  of 
rotation  of  the  sensor,  and,  in  the  last  3  columns,  provides 
the  ordered  Euler  angle  components  that  form  the  single 
maximum  rotation  angle.  The  second  file,  DISPLACE.DAT,  is 
also  an  N  X  4  matrix.  The  first  column  is  the  magnitude  of 
displacement  of  the  array,  while  the  last  3  columns  give  the 
X,  Y,  and  Z  axis  components  of  displacement. 


V.  INTERFftCE  WITH  EXISTING  PROGRftMS 

The  so-ftMare  developed  to  compute  residuals  and  generate 
simulated  track  segments  are  the  original  work  o-f  the 
author,  aided  by  such  canned  routines  as  eigensystem 
analysis,  random  number  generators  and  vector  arithmetic. 
These  canned  subroutines  came  -from  the  IMSL  Library  o-f 
Nathematic  and  Statistical  -functions  developed  for  the  IBM 
PC  computers  CRe-f.43,  and  are  compiled  in  machine  language 
libraries  not  reproducible  here.  The  author's  FORTRAN  code 
is  listed  as  Appendix  A. 

In  order  to  produce  the  estimates  o-f  sensor  displacement 
and  rotation,  however,  it  was  necessary  to  inter-face  with  a 
large  stand  alone  FORTRAN  package,  KEYMAIN.  This  program 
takes  as  input  crossover  region  data  sets  <real  or 
simulated)  and  produces  the  estimates  o-f  displacement  and 
rotation  o-f  the  second  sensor  o-f  the  crossover  pair,  based 
on  the  assumed  accurate  position  o-f  the  -first.  The 
orientation  correction  output  by  KEYMAIN  is  actually  a  three 
valued  vector  o-f  ordered  angular  rotations  in  the  XY,  XZ  and 
YZ  planes.  Similarly,  the  location  correction  is  a  three 
valued  vector  o-f  displacements  in  the  X,  Y  and  Z  directions. 
To  reduce  complexity,  the  ordered  Euler  angles  were  reduced 
to  a  single  maximum  angle  o-f  rotation  about  an  appropriately 
tilted  axis,   and  the  three  components  of  displacement   were 
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reduced  to  a  single  quantity,  magnitude  o-f  displacement. 
Thus  the  six  dimensional  quality  o-f  position  corrections  Mas 
reduced  to  tMO. 

KEYMAIN  Mas  designed  as  a  stand  alone  product  and  as 
such  is  very  user  -friendly.  It  Mas  also  designed  to  take 
not  just  one,  but  several,  crossover  data  sets  and  produce 
displacement  and  rotation  estimates  -for  several  sensors  at 
once.  Because  its  use  in  the  simulation  Mas  to  process  a 
single  simulated  crossover  data  set,  and  because  it  was 
being  called  as  a  subroutine  rather  than  used  as  a  stand 
alone  package,  signi-ficant  changes  Mere  required  in 
the  program  to  obtain  -fast  simulated  results  uninterrupted 
by  the  noM  unnecessary  user  -friendliness.  Not  only  did  this 
require  the  modi-f ication  o-f  the  executive  driver  routine, 
but  also  modi -f  ication  o-f  several  o-f  the  called  subroutines. 
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VI.  RESULTS 

It  is  imperative  that  the  simulated  track  segments 
"look"  and  "act"  like  the  original  track  segment  they  are 
simulating.  That  the  simulated  track  segments  act  like  the 
original  is  guaranteed  by  the  equations:  the  residuals  will 
have  mean  zero  by  definition  and  their  covariance  matrices 
were  computed  to  be  the  same  as  the  original's.  For  visual 
verification,  Appendix  D  is  included.  Appendix  D  contains 
plots  o-f  each  original  track  segment  overlaid  by  one  o-f  the 
tracks  simulated  from  it.  A  good  fit  is  apparent  in  all 
cases. 

The  simulation  model  can  produce  the  data  needed  to 
produce  a  scatterplot  of  magnitude  of  displacement  (in  feet) 
versus  maximum  angle  of  rotation  (in  radians)  like  the 
schematic  shown  in  Figure  5.  This  represents  an  imaginary 
situation  where  a  target  vehicle  was  driven  through  a 
single  crossover  region  on  a  range  700  times  over  the  exact 
same  track,  and  the  results  were  fed  into  the  program  to 
produce  the  displacement  and  rotation  values.  It  acts  as  a 
data  base,  or  sampling  distribution,  for  an  array  whose 
position  has  not  changed,  and  the  graph  depicts  the  natural 
variability  of  the  data.  The  contours  represent  some 
theoretical  confidence  levels  for  that  particular  sensor. 
For  example,   take  the  outermost  contour,  labeled  .90.   In  a 
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-future  tracking  exercise,  crossover  data  can  be  -fed  into  the 
program  that  produces  displacement  and  rotation  estimates. 
That  Mill  produce  a  point  on  the  graph.  If  that  point  lies 
outside  the  contour,  Me  have  good  evidence  that  the  sensor 
has  moved.  Statistically  speaking,  one  Mould  have  a  10 
percent  chance  of  rejecting  a  true  null  hypothesis,  given  a 
null  hypothesis  of  no  sensor  movement.  Conversely,  if  the 
point  plotted  closer  to  the  middle  of  the  graph,  there  Mould 
be  insufficient  evidence  to  support  the  hypothesis  of  sensor 
movement,  and  the  apparent  movement  Mould  be  attributed  to 
the  inherent  variability  in  the  data. 

It  is  unreasonable  to  expect  that  a  target  vehicle  could 
be  driven  along  the  same  track  700  times,  nor  Mould  the 
range  operators  be  likely  to  attempt  it.  The  simulation 
program,  hoMever ,  needs  only  one  straight  line  segment  of 
track  through  a  crossover  region,  and  can  then  generate  as 
many  track  data  sets  as  needed  to  produce  the  graph. 

Figures  6  through  11  are  plots  produced  from  the 
simulation  model  using  the  data  sets  from  Figure  3  (p.  18). 
Figure  9  is  a  "Mel  1 -behaved"  plot  that  could  conceivably, 
Mith  many  more  runs,  yield  the  type  of  graph  displayed  in 
Figure  5.  In  fact,  the  points  appear  so  evenly  distributed 
that  one  could  conceivably  assume  independence  betMeen  the 
rotation  and  displacement  values. 

Although  computationally  attractive,  independence  is 
thought   to  be  a  poor  assumption,   because  if  a   sensor   has 
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displaced  a  great  distance,  it  has  also  had  ample 
opportunity  to  rotate.  This  observation  is  speculative,  but 
seems  to  be  borne  out  in  the  remaining  figures  where  a 
strong  positive  correlation  appears  to  exist  between  the 
displacement  and  rotation  values.  The  fishhook  appearance 
on  these  graphs  is  a  result  of  making  the  correction 
estimates  positive  regardless  of  the  direction  of 
displacement  or  rotation.  The  graphs  taken  as  a  whole  serve 
to  illustrate  the  variable  nature  of  the  track  data. 

These  three  graphs,  drawn  from  the  simulation,  vividly 
illustrate  the  need  for  further  study.  Two  points 
concerning  the  data  on  which  these  graphs  were  based  sre 
perhaps  of  some  importance  and  may  begin  to  explain  the 
strange  nature  of  the  graphs  produced.  First,  the  data  is 
fairly  old,  taken  in  September,  1982.  The  range  operators 
from  NUWES  state  that  their  capabilities  have  improved  in 
the  interim  3  years  and  that  newer  data  could  prove 
significantly  more  accurate.  Second,  the  data  was  taken  on 
a  single  day,  from  a  single  range,  using  only  3  of  the  24 
sensors  on  the  range. 

The  great  variety  of  graphs  produced  points  to  the 
necessity  to  do  much  more  research  and  to  the  utility  of  a 
simulation  model  to  accomplish  it.  In  addition  to  the 
aforementioned  use  of  newer  data,  the  following 
considerations  warrant  study  to  ascertain  their  possible 
effects: 
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(a)  Day  to  day  variations  on  the  range  —  It  is 
not  clear  at  this  time  whether  the  character  o-f 
the  water  column  in  which  the  target  operates 
is  invariant  over  time.  Di-f -f erences  in  salinity 
and/or  temperature  could  conceivably  evolve  over 
time,  a-ffecting  the  accuracy  o-f  the  track  data. 

(b)  Seasonal  variations  on  the  range  —  Mhile  it 
is  not  clear  whether  changes  occur  on  the  range  on 
a  daily  basis,  it  certainly  seems  reasonable  to 
expect  a  variation  -From  season  to  season.  Our 
data,  taken  -from  a  single  range  on  a  single  day, 
was  insu-f -f icient  to  explore  this. 

<c>  Location  on  the  range —  Current  practice  on  a 
tracking  day  is  to  take  one  sample  o-f  the  water 
column  on  the  range,  checking  -for  salinity, 
temperature,  and  other  -factors  that  a-f-fect  the 
sound  velocity  pro-file,  and  assume  that  the 
results  hold  true  -for  the  duration  o-f  the  exercise 
•for  the  entirety  o-f  the  range.  The  possibility  o-f  a 
daily  change  was  discussed  earlier.  Here  we  mention 
the  possibility  o-f  di-f-ferent  water  characteristics 
■from  one  end  o-f  the  range  to  the  other. 

<d)  Geometry  of  tracking  runs  —  It  is  possible 
that  varying  the  geometry  of  the  tracking  run 
could  result  in  changes  to  the  quality  of  the 
data  recorded.  It  can  be  seen  from  figure  3  (on 
page  18)  that  all  of  the  data  used  for  the 
simulations  was  from  tracks  that  run  predominantly 
downrange  with  relatively  little  crossrange  change 
and  virtually  no  depth  change.  (Although  depth 
cannot  be  seen  from  figure  3,  it  was  examined  for  all 
the  tracks. ) 

<e)   Depth   of  target   —  Depth  appears  to  be   the   least 

reliable  of   the  three   dimensions   recorded   during 

tracking  runs.     Our   data   is   from   targets   that 

operated  in   a   narrow   depth   band.     Deeper    or 

shallower  targets   could   significantly   affect   the 
data. 

<f>  Inhomogenei ties  in  the  water  —  It  is  quite 
conceivable  that  undetected  inhomogenei ties  in  the 
water  could  significantly  affect  the  quality  of 
the  data.  It  would  appear  that  currents  and 
turbulence  could  affect  the  passage  of  sound 
through  the  water  column,  but  quantifying  these 
disturbances  could  be  a  major  practical  problem. 


Precisely  characterizing  the  variability  o-f  the 
correction  estimates  is  elusive  and  the  preceding 
considerations  invite  much  -future  research  be-fore  that  goal 
is  reached.  The  existence  o-f  this  simulation  model 
brightens  the  outlook  for  ultimate  success. 
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VII.  AREAS  FOR  FUTURE  WORK 

Problems  encountered  in  the  development  o-f  the 
simulation  model  and  intuition  gained  as  a  result  o-f  working 
on  it  gave  birth  to  several  areas  -for  potential  -future 
study.   Some  o-f  these  are  listed  below. 

A.  RESIDUALS  FDR  fcURVED  CROSSOVER  DATA  TRACKS 

It  Mas  convenient  in  this  simulation  model  to  use 
straight  line  segments  o-f  track  to  get  residuals.  This  is 
because  the  method  o-f  principle  components  Morks  only  -for 
straight  lines  in  N-space,  rather  than  -for  curves.  I^  a 
method  could  be  developed  to  regress  the  data  -for  any  track 
onto  its  best  path,  regardless  o-f  curvature,  one  major 
obstacle  in  the  use  of  the  simulation  Mould  be  removed. 
Whereas  noM  we  are  limited  in  the  type  o-f  data  we  can  use, 
such  a  method  would  enable  the  use  o-f  all  crossover  data. 

B.  WEIGHTED  DATA  BASED  ON  DISTANCE  TO  SENSOR  ARRAY 

It  was  assumed  in  the  model  that  recorded  track  data 
was  uni-formly  accurate  and  reliable,  regardless  o-f  the 
distance  -from  target  to  sensor.  That  is,  no  distinction  was 
made  between  the  accuracy  o-f  the  data  when  a  target  vehicle 
was  "close"  to  a  sensor  and  when  the  target  was  -far  away. 
Some  sort  o-f  weighting  scheme  may  prove  beneficial  and  help 
smooth  out  the  current  data  discrepancies. 
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C.  TR I SENSOR  CROSSOVER  DATA 

Although  smaller  in  size,  there  exist  areas  on  the 
range  where  a  target  may  be  tracked  simultaneously  by  three 
sensors  at  once.  These  tri sensor  crossover  regions  could 
provide  insight  into  the  accuracy  of  the  position  o-f  a 
sensor  and  may  yield  more  accurate  estimates  o-f  displacement 
and  rotation. 

D.  TRACK  DATA  SELECTION  PROCEDURES 

It  is  not  now  known  whether  a  greater  number  o-f  data 
points  in  a  track  segment  yields  better  results  in  the 
simulation  and  intuition  is  limited  on  this  point.  Research 
in  this  area  could  provide  valuable  guidelines  -For  selection 
o-f  data  to  use  in  the  simulation  in  the  -future. 

E.  DECISION  RULES  TO  DETERMINE  SENSOR  MOVEMENT 

A-fter  accurate  characterization  o-f  the  variability  of 
track  data  and  correction  estimates  is  made,  the  next 
logical  and  extremely  useful  step  would  be  the  construction 
of  statistically  sound  decision  rules  to  determine  the 
f ol lowing: 

(a)  The  discrepancy  noted  between  the  tracks  in  a 
crossover  data  set  can  be  explained  by  the  inherent 
variability  of  the  data. 

(b)  The  discrepancy  noted  between  the  tracks  in  a 
crossover  data  set  can  be  quantified  by  the  sensor 
slippage  model  and  can  be  computationally  corrected. 

(c)  The  discrepancy  noted  between  the  tracks  in  a 
crossover  data  set  cannot  be  explained  by  the  model 
of  sensor  movement  and  some  other  explanation  must  be 
sought. 


Clearly,  the  last  option  is  the  least  desirable,  but  if 
that  situation  is  present,  it  needs  to  be  noted. 

F.   COMPUTATIONAL  RANGE  RESURVEY 

I-F  the  method  o-f  providing  correction  estimates  can  be 
proven  to  be  reliable  and  accurate,  it  provides  a  potential 
method  o-f  range  resurvey  that  has  several  advantages  over 
the  current  method.  First,  the  range  Mould  not  need  to  be 
shut  doMn  to  resurvey.  In  -fact,  range  use  Mould  be 
mandatory  to  keep  up  to  date  sensor  array  positions. 
Second,  it  Mould  be  less  expensive  than  the  current  method, 
replacing  the  equipment  and  manpower  intensive  current 
process  Mith  relatively  inexpensive  computer  assets.  Third, 
it  Mould  take  less  time.  Resurvey  o-f  a  single  sensor  array 
on  the  range  can  take  up  to  a  day;  generating  corrections 
•from  track  data  takes  seconds.  Fourth,  since  the  current 
method  uses  a  cra-ft  on  the  sur-face  equipped  Mith  a  pinger, 
all  the  pings  must  travel  through  the  -first  150-200  -feet  o-f 
the  Mater  column,  Mhere  the  sound  velocity  pro-file  is  quite 
variable  and  most  di-fficult  to  determine,  to  get  to  the 
sensor  array.  In  contrast,  the  underMater  target  vehicles 
tracked  by  the  arrays  are  typically  in  400-600  -feet  of  Mater 
Mhere  the  sound  velocity  pro-file  is  much  smoother  and  easier 
to  predict.  Thus,  one  source  o-f  variability  in  determining 
sensor  position  is  reduced. 

It  is  not  envisioned  that  this  computer  method  could 
ever  replace  the  current  survey  process,   but  rather  augment 
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it.  Some  way  to  determine  the  position  o-f  one  sensor  is 
necessary  be-fore  KEYMAIN  can  even  begin  to  -function. 
However,  i-f  this  computer  process  could  augment  current 
resurvey  e-f-forts,  or  reduce  the  -frequency  with  which 
resurveys  must  be  made,  significant  savings  in  time  and 
money  could  result. 


APPENDIX  A 


FORTRAN  LISTING  FOR  PROGRAM  SIMDAT 


PROGRAM  SIMDAT2 
C 

C...Thi5  program  simulates  3— D  track  data  based  on  a 
C. .. specif ic  real  track  specified  by  the  user. 
C...User  inputs  are; 

C.  .  .       1.   track  segment  data  -file  to  be  simulated 
C. . .      2.   number  of  the  left  and  right  sensor  in  the 
C. . .  crossover  pair 

C. . .      3.   number  of  simulated  tracks  desired 
C. . .      4.   request  for  sample  simulated  track  (YES  or  NO) 
C. . .      5.   random  number  generator  seed 
C...The  program  output  is: 

C. . .      1.   file  of  residuals  from  the  original  track 
C...  (RESIDUAL.DAT) 

C,  .  .      2.   file  of  a  simulated  crossover  track,  if 
C...  requested  (SIMTRACK.DAT) 

C.  .  .      3.   file  of  displacement  values  (in  feet)  for  the 
C. . .  right  sensor  of  each  simulated  track  in  4 

C. . .  columns 

C. . .  Col  1     magnitude  of  displacement 

C. . .  Col  2-4   X,Y,Z  components  of  displacement 

C.  .  .      4.   file  of  rotation  values  (in  radians)  for  the 
C. . .  right  sensor  of  each  simulated  track  in  4 

C. . .  columns 

C. . .  Col  1     maximum  angle  of  rotation 

C. . .  Col  2—4   ordered  Euler  angles  of  rotation 

C 

C. . .VARIABLE  DECLARATION 
C 

SIMS, 


INTEGER*4  N,  I,  J,  K, 

lER, 

BI61 (3) ,  BI62(3) , 

POINT (130),  TRACKS, 

IDL, 

IDR,  TRKOUT 

CHARACTER  DSNAME*13 

REAL*4  NORM (260) 

REAL*8  TRACK(130,6) ,  TRl (130,3),  TR2(130,3),  MBAR1(3), 

+  MBAR2(3),  MBARMl (130,3) ,  MBARM2 ( 130 , 3) ,  COVl (3,3), 

+  C0V2(3,3),  DIFl,  PI (3,3),  MUTl (130,3),  SIM1(3,130), 

+  DDl (3) ,  DD2(3),  PA(3,3),  PB(3,3),  W0RK(130),  Dl,  D2, 

+  DIF2,  SUMl,  SUM2,  Al ,  A2 ,  Zl(3,130),  Z2(3,130), 

+  P2(3,3),  ZTl (130,3),  ZT2(130,3),  TBAR,  ZIBAR,  Z28AR, 

+  rRKSUM(130),  TKSUM2,  CT1(3,130),  CT2(3,130),  SEED, 

-t-  MUT2(130,3),  RESIDl  (130,3)  ,  RESID2  ( 130, 3)  ,  Sy2(3,3), 
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c 

c. 

c 


+   C0V1R<3,3),  CaV2R<3,3),  301(3,3),  SIMTRK (20U ,6) , 

■•-   SIM2(3,130),  ROTATE  ( 1000,4)  ,  DISP  ( 1000  ,  4)  ,  DAIA(2,4) 

Begin  the  user  input  section 


C 
C.  , 

C 


10 

c 

w  ■  a  I 

C 


20 
30 
C 

w  •  •  • 

\m^  m     «  I 
w  a  ■  ■ 

C 


WRITE(*,* 

WRITE<*.* 

READ  <*, ' 

WRITE<*,* 

WRITE(»,* 

WRITE<».* 

READ(*,*) 

WRITEC*,* 

WRITEC*,* 

READ<*,*) 

WRITE(*,* 

WRITE(*,* 

WRITE<*.* 

READ<*,*) 

WRITEC*,* 

WRITE(*,* 

WRI rE<*,* 

READ<*,*) 
WRITE <*,* 
WRITE(*,* 
WRITE<*,* 
READ(*,*)  SEED 
WRITE<*,*)  '  • 


'Enter  FNAME.FT  of  crossover  data  set 
on  disk  : 
A) • )  DSNAME 


'What  is  the  NUMBER  o-f  the  le-ft  sensor  in 
the  crossover  pair  ? 
IDL 

'What  is  the  NUMBER  o-f  the  right  sensor  ? 
IDR 

'How  many  simulated  tracks  do  you  desire 
' (NOTE  :  max  1000)  ' 
TRACKS 

'Do  you  want  a  sample  simulated  track  ?' 
'Enter  1  for  YES,  O  (zero)  for  NO  ' 
TRKOUT 

'Seed  for  the  random  number  generator' 
'NOTE  :  Seed  must  include  a  decimal 


Read  data  from  file 

OPEN  ( 1 , FILE=DSNAME , STATUS= ' OLD ' ) 

N  =  O 

N  =  N  +  1 

READ(1,*,END=30,  ERR=30)  POINT (N) , (TRACK (N, I) , 1=1 ,6) 

Separate  into  "left"  and  "right"  sensor  tracks 

DO  20  I  =  1,3 

TRl (N, I)  =  TRACK (N, I) 

TR2(N,I)  =  TRACK(N,I+3) 
CONTINUE 
GOTO  10 

CLOSE  (UNIT  =  1) 
N  =  N  -  1 

Compute  the  covariance  matrix  for  each  track. 
-First  step,  get  column  averages  (with  first  column 
average,  TBAR,  computed  for  later  use) 

TBAR  =  O. 

DO  50  I  =  1,3 
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40 


50 


C 

C. 

c. 
c. 
c. 
c 


60 


70 
80 


MBARl <I)  =0. 

MBAR2<I)  =  O. 

DO  40  J  =  1,N 

MBARl (I)  =  MBARl (I) 
MBAR2<I)  =  MBAR2(I) 


+  TRl (J, I) 
+  TR2<J,I) 


IF  (I 
CONTINUE 
MBARl <I) 
MBAR2  < I ) 

CONTINUE 

TBAR  =  TBAR 


.EQ.  1)  TBAR  =  TBAR  +  DBLE (POINT ( J ) ) 


=  MBARl (I) 
=  MBAR2<I) 

/  DBLE(N> 


DBLE(N) 
DBLE<N) 


Do  the  matrix  multiplication  :  A(t)xA  , 
the  mean  -from  each  column  entry  to  form 


subtracting  off 
the  covariance 


matrix.   Also  make  the  matrix  o-f 
later  use. 


(points  -  means)  -for 


90 

100 
C 
C. 


DO  80  I  =  1,3 
DO  70  J  =  1 , 3 
COVl (I, J)  =  O. 
C0V2(I,J)  =  O. 
DO  60  K  =  1,N 
C0V1(I,J)  =  COVl  (I,J)-KTR1  (K,I)-MBAR1  (I)  ) 
t  *(TR1 (K,J)-MBAR1 (J) ) 

CaV2(I,J)  =  C0V2(I,J)+(TR2(K,I)-MBAR2(I) ) 
ir  *(TR2(K,J)-MBAR2(J)  ) 

CONTINUE 
COVl (I, J) 
C0V2(I,J) 
CONTINUE 
CONTINUE 
DO  1 00  I  =  1 , N 
DO  90  J  =  1 , 3 
MBARMl (I, J) 
MBARM2(I,J) 
CONTINUE 
CONTINUE 


=  COVl (I ,J) 
=  C0V2(I,J) 


DBLE(N-l) 
DBLE(N-l) 


=  TRl (I, J) 
=  TR2(I,J) 


-  MBARl (J) 

-  MBAR2(J) 


C. 

C. 

C. 

C. 

C. 

C 

C. 

C. 

C 


C 
C, 


.Form  the  matrix  P  of  ordered  principle  components  for 
.each  track.   The  columns  of  P  are  the  eigenvectors 
.associated  with  the  eigen-values  of  the  covariance 
.matrix  for  each  track  arranged  in  order  of  descending 
.eigenvalues.   (ie.  the  eigenvector  associated  with 
.the  largest  eigenvalue  is  the  first  column) 

.Call  routine  to  compute  eigenvalues/vectors  for  each 
. track 

CALL  EIGRS(C0V1 ,3, 11 ,DD1 ,PA,3, WORK, lER) 
CALL  EIGRS(C0V2,3, 11 , DD2 , PB , 3 , WORK, lER) 

.Get  eigenvectors  in  eigenvalue  order,  largest  to 
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C. . . smal lest ,  by  column 

C 

C...DD1/2  =  eigenvalue  vector 

C...PA/PB  =  eigenvector  matrices 

C 

BIGl(l)  =  1 
BIG2(1)  =  1 
BIGl (3)  =  3 
BIG2<3)  =»  3 

IF  (DDKBIQKD)  .LT.  DDl  (2)  )  BIGKl)  =  2 
IF  <DD2(BIG2<1))  .LT.  DD2(2))  BIG2<1)  =  2 
IF  (DDl  (BIGl  <1)  )  .LT.  DDl  (3)  )  BIGKl)  =  3 
IF  (DD2(BIG2(1) )  .LT.  DD2(3))  BIG2(1)  =  3 
IF  (DDl (BIGl (3) )  .GT.  DDl(l))  BIGl (3)  =  1 
IF  (DD2(BIG2(3) >  .GT.  DD2 ( 1 ) )  BIG2(3)  =  1 
IF  (DDl (BIGl (3) >  .GT.  DDl (2) )  BIGl (3)  »  2 
IF  (DD2(BIG2(3> )  .GT.  DD2(2>)  BIG2(3)  =  2 
IF  ((BIGKl)  +  BIG1(3))  .EQ.  3)  THEN 

BIGl (2)  =  3 
ELSE 

IF  ((BIGKl)  +  BI61(3)>  .EQ.  4)  THEN 

BIGl (2)  =  2 
ELSE 

BIGl (2)  =  1 
END  IF 
END  IF 
IF  ((BIG2(1)  +  BIG2(3))  .ED.  3)  THEN 

BIG2(2)  =  3 
ELSE 

IF  ((BIG2(1)  -•-  BIG2(3))  .EQ.  4)  THEN 

BIG2(2)  =  2 
ELSE 

BIG2(2)  =  1 
END  IF 
END  IF 
DO  130  I  =  1,3 

DO  120  J  =  1,3 

PI  (I, J)  =  PA(I,BIG1 (J) ) 
P2(I,J)  =  PB(I,BIG2(J) ) 
120      CONTINUE 
130   CONTINUE 
C 

C... Compute  the  matrix  ZT  ior    each  track 

C...ZT  represents  the  projection  of  the  track  data  onto 
C...the  principle  components 

C ZT  =  (TR  -  MBAR)  x  P  =  MBARM  x  P 

C... where  TR  -  MBAR  =  track  data  minus  the  column  average 

C...for  each  row 

C 

C...Call  routine  to  multiply  matrices  AxB 

C 

CALL  VMULFF(MBARM1 ,P1 ,N,3,3, 130,3,ZT1 , 130, lER) 
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CALL  VMULFF ( MBARM2 , P2 , N , 3 , 3 , 1 30 , 3 , ZT2 , 1 30 , I ER ) 
C 

C... Since  there  are  some  points  missing  -from  the  data  set, 
C... per form  a  simple  least  squares  linear  regression  onto 
C...the  -first  principle  component 
C 

TKSUM2  =0.0 
ZIBAR  =  0.0 
Z2BAR  =0.0 
DO  140  I  =  1,N 

ZIBAR  =  ZIBAR  -»•  ZT1<I,1) 
Z2BAR  =  Z2BAR  ■•■  ZT2  <  1 , 1 ) 
TRKSUM(I)  =  DBLE(POINT<I) )  -  TBAR 
TKSUM2  =  TKSUM2  +  TRKSUM < I > **2 
140   CONTINUE 

ZIBAR  =  ZIBAR  /  DBLE(N) 

Z2BAR  =  Z2BAR  /  DBLE<N) 

SUMl  =0.0 

SUM2  =0.0 

DO  150  I  =  1,N 

DIFl  =  <ZT1(I,1)  -  ZIBAR)  »  TRKSUM < I) 
DIF2  =  <ZT2<I,1)  -  Z2BAR)  *  TRKSUM < I) 
SUMl  =  SUMl  +  DIFl 
SUM2  =  SUM2  +  DIF2 
150   CONTINUE 

Dl  =  SUMl  /  TKSUM2 
D2  =  SUM2  /  TKSUM2 
Al  =  ZIBAR  -  Dl  *  TBAR 
A2  =  Z2BAR  -  D2  *  TBAR 
DO  160  I  =  1,N 

ZT1(I,1)  =  Al  +  Dl  *  DBLE(POINT<I) ) 
ZT2<I,1)  =  A2  +  D2  *  DBLE<POINT<I) ) 
160   CONTINUE 
C 

C...Get  CT  matrix  which  represents  the  orthogonal  projection 
C...of  the  data  onto  the  straight  line  of  the  first 
C. .. principle  component 
C 

DO  170  I  =  1,N 

Zl  <1,I)  =  ZTl (1,1) 
Z2(1,I)  =  ZT2(I,1) 
170   CONTINUE 

DO  ISO  I  =  1,N 
Zl (2,1)  =  0.0 
Zl (3,1)  =  0.0 
Z2(2,I)  =  0.0 
Z2(3,I)  =  0.0 
180  CONTINUE 
C 

C...Call  routine  to  multiply  matrices  AxB 
C 

CALL  VMULFF(P1,Z1,3,3,N,3,3,CT1,3,IER) 
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c 

c, 

c 


CALL  VMULFF (P2,Z2,3,3,N,3,3, CT2 , 3 , I ER ) 
Move  "line"  o-f  data  back  into  original  coordinats  system 


DO  200  I  =  1,N 

DO  190  J  =  1,3 

MUTl (I, J)  =  CTl (J, I) 
MUT2<I,J)  =  CT2<J,I) 


MBARl (J) 
MBAR2(J) 


C... Compute  residuals  for  each  track 
C 

RESIDl (I, J)  =  TRl (I, J) 
RESID2<I,J)  =  TR2<I,J) 
190      CONTINUE 
200   CONTINUE 
C 

Write  the  set  o-f  residuals  out 


-  MUTl (I, J) 

-  MUT2<I,J) 


C. 
C 


210 
C 
C. 


to  the  -file  RESIDUAL.DAT 


OPEN  (2,  FILE  =  'RESIDUAL.DAT',  STATUS  =  'NEW') 
DO  210  I  =  1,N 

WRITE (2,330)  (RESIDl ( I , J) , J=l ,3) , (RESID2(I,J) ,J=1 ,3) 
CONTINUE 


C. 

C 

C. 

C. 

C 


.Compute  the  covariance  matrix  of  the  residuals 
.(Note  :  column  averages  are  identically  zero) 


.Call  routine  to  multiply  matrices  A(t)xB, 
.by  N-1 


then  divide 


CALL  VMULFM(RESID1,RESID1,N,3,3,130, 130,C0V1R,3, lER) 
CALL  VMULFM(RESID2,RESID2,N,3,3, 130, 130,CGV2R,3, lER) 
DO  230  1=1,3 

DO  220  J  =  1 ,3 

=  C0V1R(I,J)  /  DBLE(N-l) 
=  C0V2R<I,J)  /  DBLE(N-l) 


C0V1R(I,J> 
CGV2R<I ,J) 
CONTINUE 
CONTINUE 


220 

230 
C 
C...Get  "square  root"  o-f  covariance  matrix  o-f  residuals 


SQl  (3 
SQ2<3 
SDl  (2 
SQ2<2 
BQl  (2 
SQ2(2 
SQl  (1 
SQ2(1 
SDl  (1 
SQ2(1 
SQl  (1 
SQ2(1 


3) 
3) 
3) 
3) 
2) 
2) 
3) 
3) 
2) 
2) 
1) 
1) 


DSQRT(C0V1R(3,3) ) 
DSQRT(CGV2R(3,3) ) 
CGV1R(2,3)  /  SQl (3,3) 
CGV2R(2,3)  /  SQ2(3,3) 
DSQRT(CGV1R(2,2)  -  SQ1(2,3)»«2) 
DSQRT(CGV2R(2,2)  -  SQ2(2,3)»*2) 
CaVlR(l,3)  /  SQl (3,3) 
C0V2R(1,3)  /  SQ2(3,3) 
(CGV1R(1,2)  -  SQl (1,3)*SQ1 (2,3) ) 
(C0V2R(1,2)  -  SQ2(1 ,3)*SQ2(2,3) ) 


SDl (2,2) 
SQ2(2,2) 


=  DSQRT(C0V1R(1,1)-(SQ1 (1,2)**2+SD1 (1,3) ♦♦2) ) 
=  DSQRT(CGV2R(1  ,  1)-(SQ2(1  ,2)  *»2-»-SD2  (  1  ,3)  ♦♦2)  ) 
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SQl <2,1)  =  0. 
502(2,1)  =  0. 
SQl (3,2)  =  O. 
SQ2(3,2)  =  0. 
SQl (3,1)  =  0. 
SQ2(3,1)  =  O. 
C 

C...  Compute  sets  o-f  residuals  and  get  rotation/displacement 
C. . . values 
C 

DO  410  SIMS  =  1, TRACKS 
C 

C.  ..Compute  set  o-f  simulated  residuals  from  normal  (0,1) 
C. . . devi  ates 
C 

DO  260  I  =  1,3 
C 

C...Call  routine  to  generate  Normal (0,1)  deviates 
C 

CALL  GGNPM( SEED, 2*N, NORM) 
DO  250  J  =  1 , N 

RESIDl (J, I)  =  NORM(J) 
RESID2(J,I)  =  NORM(N+J) 
250      CONTINUE 
260   CONTINUE 
C 

C...Call  routine  to  multiply  matrices  AxB(t) 
C 

CALL  VMULFP(SQ1,RESID1,3,3,N,3,130,SIM1 ,3,IER) 
CALL  VMULFP<SQ2,RESID2,3,3,N,3,130,SIM2,3,IER) 
C 

C...Put  together  Nx6  matrix  of  simulated  tracks  for  both 
C... arrays  by  adding  straight  line  in  original  coordinate 
C... system  to  residuals 
C 

DO  280  I  =  1,N 

DO  270  J  =  1,3 

SIMTRK(I,J)  =  SIM1(J,I)  +  MUT1(I,J) 
SIMTRK(I,J+3)  =  SIM2(J,I)  +  MUT2(I,J) 
270      CONTINUE 
280   CONTINUE 
C 

C...l«Jrite  the  first  simulated  track  out  to  the  file 
C. . . SIMTRACK. DAT  if  a  sample  simulated  track  was  requested. 
C 

IF ((SIMS  .EQ.  1)  .AND.  (TRKOUT  .GT.  O) )  THEN 

OPEN  (3,  FILE  =  'SIMTRACK.DAT',  STATUS  =  'NEW') 
DO  285  I  =  1,N 

WRITE(3,340)  (SIMTRK ( I , J) , J  =  1,6) 
285      CONTINUE 
END  IF 
C 
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C...Feed  the  simulated  track  into  KEYMAIN  to  get  rotation 

C...and  displacement  numbers 

C 

CALL  KE YSUB ( S I MTRK , N , DATA , I DL , I DR ) 
C 

C ..Make  2  matrices  -  one  for  displacement  data  and  one  for 
C...the  rotation  data 
C 

DO  290  I  =  1,4 

DISP(SIMS,I)  »  DATA < 1,1) 
ROTATE (SIMS, I)  =  DATA (2, I) 
290   CONTINUE 

WRITE(*,320)  SIMS 
C 

C...Go  back  and  do  it  again 
C 

410   CONTINUE 
C 

C. ..After  TRACKS  simulated  tracks,  write  the  displacement 
C...set5  and  the  rotation  sets  out  to  a  file 
C 

OPEN (4, FILE  =  'DISPLACE.DAT' , STATUS  =  ' NEW ' ) 
OPEN <5, FILE  =  'ROTATE.DAT' , STATUS  =  'NEW') 
DO  300  1=1, TRACKS 

WRITE(4,310)  <DISP<I,J),J  =  1,4) 
WRITE(5,310)  <ROTATE<I,J) ,J  =  1,4) 
300   CONTINUE 
C 

C... Close  out  the  files 
C 

CLOSE (UNIT  =  2) 
CLOSE (UNIT  =  3) 
CLOSE (UNIT  =  4) 
CLOSE (UNIT  =  5) 
STOP 
310   F0RMAT(2X,4F17.a) 

320   FORMAT (2X, 'Through  KEYMAIN  ',16,'  time(s)  so  f ar ' ) 
330   F0RMAT(1X,6F12.7) 
340   F0RMAT(1X,6F12.5) 
END 
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APPENDIX  E 
FORTRAN  SUBROUTINES  CONECT  AND  REDUCE 

The  program  KEYMAIN  requires  that  any  array  for  which 
correction  estimates  are  desired  must  be  "connected"  to  the 
first  input  sensor  array  by  one,  or  a  series  of,  crossover 
data  sets.  For  example,  if  an  input  crossover  data  set  uses 
sensor  arrays  5  and  6,  while  another  uses  arrays  6  and  10, 
all  three  of  the  arrays  <5,  6  and  10)  are  connected.  If, 
however,  a  first  input  data  set  uses  arrays  7  and  9,  a 
second  input  data  set  uses  12  and  10,  while  a  third  input 
data  set  uses  9  and  13,  the  arrays  7,  9  and  13  are 
connected,  12  and  10  are  connected,  but  all  five  arrays  do 
not  form  a  single  connected  set.  In  this  case,  if  7  was  the 
first  array  input  to  the  program,  correction  estimates  could 
not  be  made  for  arrays  10  and  12.  The  subroutine  CONEC I 
checks  to  see  that  connectedness  exists  in  the  input  data 
before  KEYMAIN  is  allowed  to  continue. 

KEYMAIN  allows  3  options  if  CONECT  discovers  that  the 
arrays  of  the  input  data  sets  Bre  not  connected.  One  option 
is  to  quit,  in  which  case  the  program  terminates.  Another 
option  is  to  add  more  data  so  that  all  the  arrays  are 
connected.  In  the  second  example  above,  for  instance,  if  a 
crossover  data  set  using  arrays  9  and  12  was  input  into  the 
program,   all   5  arrays  would  then  be  connected  and   KEYMAIN 
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Mould  continue.  A  third  option  is  to  continue  with  KEYMAIN, 
but  use  only  the  first  connected  set,  that  is,  all  the 
arrays  that  are  connected  to  the  le-ft  array  o-f  the  -first 
input  crossover  data  set.  This  option  presents  special 
problems  because  it  requires  a  reduction  o-f  the  data 
structures  that  have  been  built  as  each  data  set  was  input. 
The  subroutine  REDUCE  does  this  data  structure  reduction  and 
is  called  from  KEYMAIN  only  when  this  third  option  is 
selected. 

KEYMAIN  passes  to  CONECT  two  pieces  of  information. 
The  first  is  the  variable  Rl  that  represents  the  number  of 
data  sets  that  were  input  into  KEYMAIN.  The  second  piece  of 
information  is  a  2  x  Rl  matrix,  IND2,  that  contains  the 
number  of  the  left  and  right  sensors  of  the  Rl  crossover 
data  sets.  Row  1  contains  the  left  sensor  numbers,  row  2 
the  right.  CONECT  performs  its  connectedness  check  by 
starting  a  variable  length  list  that  contains  the  array 
numbers  of  those  arrays  that  are  connected.  The  list  starts 
with  only  two  entries,  the  left  and  right  sensor  arrays  of 
the  first  crossover  data  set.  These  are  connected,  and  the 
left  sensor  is  the  "root"  to  which  all  arrays  should 
connect.  Elements  are  added  to  the  list  by  sequencing 
through  the  list  from  the  beginning,  and  adding  to  the  list 
any  array  numbers  for  IND2  that  (1)  are  not  yet  on  the  list 
and  (2)  are  connected  by  a  crossover  data  set  to  an  array 
that  is  already  on  the  list.    If,   after  sequencing  through 
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the  variable  length  list,  there  are  arrays  remaining  in  IND2 
that  never  were  put  on  the  list,  those  arrays  were  not  a 
part  o-f  the  first  connected  set.  If  all  arrays  in  IND2  were 
included  on  the  list,  all  the  arrays  were  connected.  After 
the  first  list  is  exhausted,  CONECT  repeats  this  procedure 
as  many  times  as  there  are  disjoint  sets  of  arrays,  starting 
each  new  list  with  the  arrays  of  a  data  set  not  in  any 
previous  set.  CONECT  informs  the  user  of  the  individual 
sets  of  connected  sets  of  arrays  and  raises  a  flag  to  alert 
the  user  if  all  sets  were  not  connected. 

If  all  the  arrays  in  the  input  crossover  data  sets  were 
not  connected,  the  user  has  three  options  to  proceed, 
described  earlier.  If  he  chooses  to  continue  using  the 
first  set  of  connected  arrays,  REDUCE  is  called  to  pare  the 
data  strauctures  built  up  during  the  data  input  process.  In 
particular,  the  Rl  x  3  x  3  array  CROSSA  and  the  Rl  x  6 
matrix  mean  need  to  be  reduced  to  contain  only  those 
elements  that  correspond  to  the  data  sets  in  the  first 
connected  set.  CROSSA,  which  has  Rl  "pages"  of  3  x  3 
matrices  of  crossproduct  deviations  from  the  mean,  needs  to 
have  those  pages  removed  that  correspond  to  every  crossover 
data  set  in  the  original  input  not  connected  to  the  first 
array.  MEAN,  which  has  Rl  rows  of  the  colunm  averages  for 
each  data  set,  needs  to  have  those  rows  removed  that 
correspond   to  data  sets  not  connected  to  the   first   array. 
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The  variable  Rl  itsel-f  will  be  reduced  to  re-Flect  this 
smaller  subset  o-f  connected  array  pairs. 

CONECT  stores  the  information  that  indicates  which  data 
sets  in  IND2  are  connected  to  the  -first  array.  This 
information  is  passed  to  REDUCE  through  KEYMAIN.  REDUCE 
re-forms  the  data  structures  to  re-flect  the  smaller  number  of 
data  sets  now  being  considered  by  KEYMAIN. 

CONECT  also  provides  the  matrix  INDl ,  a  K  x  Rl  matrix 
that  is  used  elsewhere  in  KEYMAIN.  The  Rl  columns  represent 
the  crossover  data  sets  input  into  KEYMAIN.  The  variable  K 
represents  the  number  of  individual  arrays  in  the  Rl  data 
sets,  each  row  representing  a  separate  array.  For  each 
column  in  INDl,  the  entries  are  all  O  except  in  the  rows 
that  represent  the  left  and  right  sensor  for  that  column's 
data  set.  If  the  row  corresponds  to  the  left  array,  the 
value  is  1.  If  it  corresponds  to  the  right  array,  the  value 
is  2.  If  REDUCE  is  called,  some  columns  (representing 
crossover  data  sets)  and  rows  (representing  individual 
sensors)  of  INDl  need  to  be  omitted.  REDUCE  reforms  INDl  to 
its  smaller  size. 
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APPENDIX  F 


FORTRAN  LISTING  FOR  SUBROUTINE  CONECT 


SUBROUTINE  CONECT  (0UT,R1 , IND2,K, INDl , IA,TESTC, IND2R, 
♦  DATSET) 

C 

C*********************»***********»**»»******************»** 
C    This  subroutine  checks  -for  the  connectedness  of  the 
C    input  data  sets.  I-f  the  problem  is  connected  then  the 
C    user  is  in-formed  and  the  array  pairs  are  printed  on  the 
C    screen;  if  not  connected,  then  the  user  is  prompted  to 
C    select  one  of  three  options  -  quit,  add  conecting  data 
C    sets,  or  run  the  program  using  the  first  connected  set 
C    that  Mas  input.  Gygax  -  July  1985 

C*«***********************»*****»****»* 

c 

C       ...Variable  declarations. 

C 

INTEGER*4  Rl ,K, IND2 (2,30) , INDl (30,30) , I , J , lA <30) , FIRST 
INTEGER*4  LIST(30) ,BEGIN,HALT ,DISCON,L,M,0,TES r C,OU ( , 
I NTEGER*4  DATSET ( 30 ) , COUNT , SAVE ( 2 , 30 ) ,  I ND2R (2,30) 

C 

C       ...Initialize  the  values  of  FIRST  and  COUNT: 

C 

FIRST  =  O 
COUNT  =  O 

C 

C       ...Make  vector  I A  =  list  of  all  arrays  (w/o  repeats) 

C  in  IND2  and  get  the  value  for  K  =  #  of  individual 

C 

C 


30 


40 
50 
60 

C 


.  , 

..Make  vector  I A  = 

=  li 

St  of 

in  IND2 

and  get 

the 

valu 

arrays. 

IA(1) 

=  IND2(1,1) 

IA(2) 

=  IND2(2,1) 

K  = 

=  3 

IF 

(R] 

L  .EQ. 

1)  GOTO 

60 

DO 

50 

I  =  1, 

,R1 

DO 

40  J  = 
M  =  K 
DO  30 

=  1,2 
-  1 
L  =  1,M 

IF 

(IND2(J, 

I) 

-EQ. 

CONTINUE 

IA(K) 

=  IND2(J 

,1) 

K  =  K 

+  1 

CONTINUE 

CONTINUE 

K  = 

=  K 

-  1 

IA(L) )  GOTO  40 
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WRITE(OUT,») 'Rl   ' ,R1 
WRITE  (OUT,*)  'K     '  ,K. 
C       ...For  each  column  o-f  INDl  (columns  correspond   to 
C  data  sets)  the  entries  are  all  zero  except  for 

C  the  row  that  corresponds  to  the  le-ft  arr^y  ''-  1) 

C  and  the  right  array  (=  2). 

C 

DO  80  I  =  1 ,R1 
DO  70  J  =  1,K 

INDl  (J,  I)  =  O 

IF  (IND2(1,I)  .EQ.  IA(J))  INDl (J, I)  =  1 
IF  (IND2(2,I)  .EQ.  IA(J))  IND1(J,I)  =  2 
70        CONTINUE 
SO    CONTINUE 
C 

C       ...Check  to  see  i-f  all  the  arrays  Are    connected. 
C 

TESTC  =  1 
LIST(l)  =  -IA(1) 
DO  131  I  =  1,R1 

IF  (IND2(1,I)  .EQ.  -LIST(l))  IND2(1,I)  =  -IND2(1,I) 
IF  (IND2(2,I)  .EQ.  -LIST(l))  IND2(2,I)  =    -IiMD2(--,I) 
131   CONTINUE 
BEGIN  =  1 
HALT  =  1 

140  IF  (.NOT.   (BEGIN  .LE.  HALT))  GOTO  170 
NODE  =  LIST(BEGIN) 

BEGIN  =  BEGIN  +  1 
DO  150  I  =  1,R1 
IF  (.NOT.  (  (NODE.EQ.  IND2(1 , I)  )  .AND.  <IND2(2, n  .GT,0^  n 

*  GOTO  150 
HALT  =  HALT  +  1 
LIST(HALT)  ^  -IND2(2,I) 
DD  141  J    =  1 , R 1 

IF  (IND2(1 ,J) .EQ.-LIST(HALT) )  IN02 ( 1 , J ) =- I ND2 ( I . J ) 
IF  (IND2(2,J)  .EO. -LIST  (HALT)  ^  IND2(2,J)  -^- IND2  (2  ,  J  ) 

141  CONTINUE 
150   CONTINUE 

OCl  1 60  I  =  1  ,  R 1 
IF  ( . NOT .  (  ( NODE . EQ .  I ND2 (2,1)). AND .  ( I ND2  ^1  , I )  . GT . O  ^  )  ? 

*  GOTO  160 
HALT  --=  HALT  +  1 

LIST (HALT)  =  -IND2(1,I) 

-L I ST ( HAL  T )  )  I ND2 (1  , J  >  =~ I ND2 ( 1 , J  > 
-LIST  (HALT)  )   IND2(2,J)-=-:ND2(2,  J) 


DO  151  J  =  1,R1 

IF  (IND2(1 ,J) .EQ 

IF  (IND2(2,J).EQ 

151 

CONTINUE 

1 60 

CONTINUE 

GOTO  140 

1  7^") 

CONTINUE 

c 

.  .  .  Pr  i  n  t  out  th^e  m 

matched  pairs, 
1  '-^  I 


DISCON  =  O 
WRITE(0UT,230) 
DO  200  I  =  1,R1 

IF  (IND2(1,I)  .LT.  O)  GOTO  190 
IF  (IND2(1,I)  .EQ.  O)  GOTO  200 

IF  ((IND2(1,I)  .GT.  0)  .AND.   (DISCON  .ED.  D)  GOTO  200 
FIRST  =  FIRST  +  1 
DISCON  =  1 
TESTC  =  0 
BEGIN  =  1 
HALT  =  1 

LIST(l)  =  IND2<1,I) 
GOTO  200 
190   WRITE(0UT,240)  -IND2 (1 , I ) , -IND2 ( 2 , I ) 

IF  (  (FIRST.  EQ.O)  .OR.  (  (FIRST.  EQ.  1)  .AND.  (DISCON.  EQ,  i)  )  ,' 
*  THEN 

COUNT  =  COUNT  +  1 

I ND2R ( 1 , COUNT )  =  - I ND2 (1,1) 

IND2R(2,C0UNT)  =  -IND2(2,I) 

DATSET (COUNT)  =  I 


END  IF 

SAVE (1,1)  =  -IND2(1,I) 

SAVE (2, I)  =  -IND2(2,I) 

IND2(1 ,1)  =  0 

IND2(2, I)  =  0 

00 

CONTINUE 

IF  (DISCON  .EQ.  1)  GOTO 

DO  220  I  =  1,R1 

IND2(1  ,1)  =  SAVEd  ,  I) 

IND2<2,I)  =  SAVE (2, I) 

20 

CONTINUE 

RETURN 

30 

FORMAT (IX,' THE  FOLLOW I NG 

140 

FORMAT (IX, 1415) 

END 

140 


PAIRS  ARE  CONNECTED 


lo: 


APPENDIX  G 


FORTRAN  LISTING  FOR  SUBROUTINE  REDUCE 


SUBROUT I NE  REDUCE  <  CROSSA , MEAN , R 1 , K , I ND 1 , I A , 
»  IND2R,DATSEI ) 

C 

c**************»***»*»******»*»*»****»»**»****»«**»»*»*»»»** 

C  This  is  a  specialized  subroutine  that  is  used   when 

C  option  three  is  invoked  as  a  result  o-F  a  -failed  con- 

C  nectedness  test.  The  disconnected  data  sets  must  be 

C  removed  from  the  variables  CROSSA  and  MEAN,  and  other 

C  program  supporting  variables  must  be  adjusted. 

C  Gygax  -  July  1985 

C 

C       ...  Variable  declarations. 

C 

INTEGER*4  Rl ,K, IND1<30,30) , lA <30) ,IND2R(2,30) ,I,J,L, 
2  M,DATSET<30) 

C 

REAL*8  CROSSA ( 30 , 3 , 3 ) , MEAN ( 30 , 6 ) 
C 

C       ...  Compute  the  new,  reduced  Rl: 
C 

DO  10  I  =  1,30 

IF  <IND2R(1,I)  .EQ.  0)  GOTO  20 
10    CONTINUE 
20    Rl  =  I  -  1 
C 

C       ...  Make  new,  reduced  vector  I A  =  list  o-f  all  arrays 
C  in  IND2R  w/o  repeats.  Also,  co(npute  a  new  K. 

C 

IA(1)  =  IND2R(1,1) 
IA<2)  =  IND2R(2,1) 
K  =  3 

IF  (Rl  .EQ.  1)  GOTO  60 
DO  50  I  =  1,R1 
DO  40  J  =  1 , 2 
M  =  K  -  1 
DO  30  L  =  1 , M 

IF  (IND2R(J,I)  .EQ.  IA(L))  GOTO  40 
30  CONTINUE 

IA(K)  =  IND2R(J,I) 
K  =  K  +  1 
40       CONTINUE 
50    CONTINUE 
60    K  =  K  -  1 
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C       ...  Remake  the  reduced  matrix  INDl  -  for  each  column 
C  in   INDl   (corressponding   to  a   data   set)   the 

C  entries  are  zero  except  for  the  entries  corres 

C  ponding  to  the  left  array  <=  1)  and  the  right 

C  array  (=  2) . 

C 

DO  80  I  =  1,R1 
DO  70  J  =  1,K 

INDl <J,I)  =  O 

IF  (IND2R(1,I)  .EQ.  IA(J))  INDl (J, I)  =  1 
IF  <IND2R<2,I)  .EQ.  IA<J))  IND1<J,I)  =  2 
70       CONTINUE 
80    CONTINUE 
C 

C       ...  Reduce  the  arrays  CROSSA  and  MEAN  to  account 
C  for  the  removed  data  sets. 

C 

DO  120  I  =  1,R1 
DO  90  J  =1,6 

MEAN (I, J)  =  MEAN(DATSET<I) ,J) 
90       CONTINUE 

DO  110  J  =  1,3 

DO  100  L  =  1,3 

CROSSA (I, J, L)  =  CROSSA<DATSET<I) ,J,L) 
100         CONTINUE 
110      CONTINUE 
120   CONTINUE 
RETURN 
END 
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